function diff = root_psi( x, ind, param)

    param.psi = x;
    x0= fzero(@ (x) (J_r_p(param.m_h+abs(x),param)), [0.1, 40], optimset('display','off'));%-(param.c+param.K)
    param.m_e=param.m_h+abs(x0);

    param.delta_2 = param.m_e.^(-1/(1-param.alpha)).*(delta_m(param.m_e,param)-(1-param.omega_1)./param.alpha./param.c.*param.m_e+(param.omega_0+param.r.*(param.c+param.K)+param.zeta.*param.c)./param.c);

    if delta_u_p(param.m_e,param)>0 || delta_m(param.m_e,param)<0 % for out of equilibrium values of psi
        param.m_u=param.m_e;
    else
        x0_u= fzero(@ (x) delta_u(abs(x)+param.m_e,param), 10^(-5), optimset('display','off'));
        param.m_u=param.m_e+abs(x0_u);
    end

    psi = q_fun(param.m_l,param)./q_fun(param.m_u,param);

    if ind == 1
        diff = (psi-param.psi);
    else
        diff = param;
    end
end